home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 5 / MacMania 5.toast / / Tools&Utilities / Plotfoil 3.2 / foil_lib.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-18  |  7.4 KB  |  294 lines  |  [TEXT/MMCC]
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <ctype.h>
  5. #include "foil_lib.h"
  6. #include "lib.h"
  7.  
  8. /*
  9.  * get_coords
  10.  */
  11. static int get_coords(FILE *fp, float *pt1, float *pt2, float *pt3)
  12. {
  13.    char buf[TXTSIZ], *p;
  14.    float x1, x2, x3;
  15.    int n;
  16.    
  17.    while (1) {
  18.       if(!fgets(buf, TXTSIZ, fp))
  19.      return EOF;
  20.  
  21.       /* skip comments */
  22.       p = buf;
  23.       while(isspace(*p)) p++;
  24.       if (in(*p, "%#\n"))
  25.      continue;
  26.         
  27.       /* scan string */
  28.       n = sscanf(buf, "%f %f %f", &x1, &x2, &x3);
  29.  
  30.       if (n == EOF)
  31.      return EOF;
  32.  
  33.       /* valid data? */
  34.       if (n > 0)
  35.      break;
  36.  
  37.       /* otherwise who knows? just keep reading */
  38.    }
  39.  
  40.    *pt1 = x1;
  41.    *pt2 = x2;
  42.    if (pt3)
  43.       *pt3 = x3;
  44.    return n;
  45. }
  46.  
  47. int read_foil(FILE *fp, char *f, airfoil_data_t *info)
  48. {
  49.    static float X[MAXPOINTS];              /* x points */
  50.    static float Y[MAXPOINTS];              /* y points */
  51.    static float YY[MAXPOINTS];             /* other y, for NACA format */
  52.    float x=0., y=0., yy=0., minx, maxx, miny, maxy;
  53.    int i, j, n, nread, npoints;
  54.    char *p;
  55.  
  56.    fgets(info->title, TXTSIZ, fp);
  57.    if((p = strchr(info->title, '\n')) != 0)
  58.       *p = 0;
  59.    n = get_coords(fp, &x, &y, &yy);
  60.  
  61.    if(n < 2 || n > 3) {
  62.       fprintf(stderr, "Syntax error in line 2 of \"%s\"!\n", f);
  63.       return 0;
  64.    }
  65.    
  66.    minx = maxx = X[0] = x;
  67.    Y[0] = y; YY[0] = yy;
  68.    i = 1;
  69.  
  70.    while(n == (nread = get_coords(fp, &x, &y, &yy))) {
  71.       X[i] = x; Y[i] = y; YY[i] = yy;
  72.       minx = min(x, minx);
  73.       maxx = max(x, maxx);
  74.             
  75.       if(i++ >= MAXPOINTS) {
  76.      fprintf(stderr,
  77.          "Too many points! Change MAXPOINTS (currently = %d) and recompile.\n", i);
  78.      return 0;
  79.       }
  80.    }
  81.  
  82.    /* any error codes? could be more robust */
  83.    switch (nread) {
  84.       /* data OK */
  85.    case 0:  break;
  86.       /* EOF */
  87.    case -1: break;
  88.       /* bad input */
  89.    case -2: fprintf(stderr, "Syntax error at \"%s\" (line %d)\n", f, i+2);
  90.       break;
  91.       /* bogus - shouldn't get here */
  92.    default:            /* NOTREACHED */
  93.       break;
  94.    }
  95.  
  96.    if(n == 3) {
  97.       /*
  98.        * NACA format; combine into one array.
  99.        */
  100.       npoints = 0;
  101.       for(j = i-1; j > npoints; j--, npoints++) {
  102.      float temp;
  103.      temp = X[npoints]; X[npoints] = X[j]; X[j] = temp;
  104.      temp = Y[npoints]; Y[npoints] = Y[j]; Y[j] = temp;
  105.       }
  106.       npoints = i;
  107.       for(j = 1; j < i; j++, npoints++) {
  108.      X[npoints] = X[i];
  109.      Y[npoints] = YY[i];
  110.       }
  111.    }
  112.    else
  113.       npoints = i;
  114.    
  115.    if(X[0] != X[npoints-1] || Y[0] != Y[npoints-1]) {
  116.       fprintf(stderr,
  117.           "Warning: data does not start and end at the same point.\n");
  118.       X[npoints] = X[0]; Y[npoints] = Y[0];
  119.       npoints++;
  120.    }
  121.    if(minx < 0.0)
  122.       for(i = 0; i < npoints; i++)
  123.      X[i] -= minx;
  124.  
  125.    if(maxx == 0.0) {
  126.       fprintf(stderr, "Error: no x values?\n");
  127.       return 0;
  128.    }
  129.  
  130.    info->npoints = npoints;
  131.  
  132.    miny = maxy = 0.;
  133.    for(i = 0; i < npoints; i++) {
  134.       info->points[i].x = X[i]/maxx;
  135.       info->points[i].y = Y[i]/maxx;
  136.       maxy = max(maxy, info->points[i].y);
  137.       miny = min(miny, info->points[i].y);
  138.    }
  139.    info->maxy = maxy;
  140.    info->miny = miny;
  141.  
  142.    return 1;
  143. }
  144.  
  145. void writefoil(FILE *fout, airfoil_data_t *foil, int tchanged, int cchanged)
  146. {
  147.    int i;
  148.    
  149.    fprintf(fout, "%s", foil->title);
  150.    if(tchanged)
  151.       fprintf(fout, " T-%.3g%%", foil->thickness * 100);
  152.    if(cchanged)
  153.       fprintf(fout, " C-%.4g%%", foil->camber * 100);
  154.    fputc('\n', fout);
  155.  
  156.    for(i = 0; i < foil->npoints; i++)
  157.       fprintf(fout, "%f %f\n", foil->points[i].x, foil->points[i].y);
  158. }
  159.  
  160. /*
  161.  * Interpolate the function at x, given 2 points on either side.  If
  162.  * straightl is set, just do a straight line (for the first intervals, when P4
  163.  * or P1 may not exist).
  164.  *
  165.  * If P1 and P4 are on opposite sides of (P2,P3) we use the straight line.
  166.  */
  167. #define lin_interp(x, x2, y2, x3, y3) ((y2)+((x)-(x2))/((x3)-(x2))*((y3)-(y2)))
  168. static double interp(double x,
  169.              double x1, double y1, double x2, double y2,
  170.              double x3, double y3, double x4, double y4, int straightl)
  171. {
  172.    double i1x, i1y, i2x, i2y;    /* the two interior points for interpolation */
  173.    
  174.    /* A is (P1,P2); B is (P2,P3); C is (P3,P4). */
  175.    point_t ma, mb, mc;
  176.    ma.x = x2 - x1; ma.y = y2 - y1;
  177.    mb.x = x3 - x2; mb.y = y3 - y2;
  178.    mc.x = x4 - x3; mc.y = y4 - y3;
  179.  
  180.    if(straightl ||
  181.       ((m_gt(ma, mb) && m_gt(mb, mc)) ||
  182.        (m_gt(mc, mb) && m_gt(mb, ma))
  183.        )) {
  184.       if((x3 - x2) < EPSILON)
  185.      /* Indeterminate */
  186.      return y3;
  187.       return lin_interp(x, x2, y2, x3, y3);
  188.    }
  189.  
  190.    i1x = x2 + (x3-x2)/3.;
  191.    i2x = x2 + (x3-x2)*2./3.;
  192.    i1y = lin_interp(i1x, x1, y1, x2, y2);
  193.    i2y = lin_interp(i2x, x3, y3, x4, y4);
  194.    if(x < i1x)
  195.       return lin_interp(x, x2, y2, i1x, i2y);
  196.    if(x > i2x)
  197.       return lin_interp(x, i2x, i2y, x3, y3);
  198.    return lin_interp(x, i1x, i1y, i2x, i2y);
  199. }
  200.  
  201. /*
  202.  * Takes an array of points, and for each each one finds the corresponding
  203.  * point on the other side by interpolating. Returns corresponding points in
  204.  * the "new" array.
  205.  *
  206.  * For now, this cheap version just finds the linear interp on the other side.
  207.  * Most airfoil data seems pretty smooth and the data points are close
  208.  * together so the linear interpolation works fine; so let's leave the "right"
  209.  * interpolation for another time. "TODO 3.?" I guess.
  210.  */
  211. void get_otherside(point_t *new, point_t *old, int n)
  212. {
  213.    int i, j, do_next, other;
  214.    /* These are for the flags: si means i is at the start; sj is j at start. */
  215.    int s, si=1, sj=1;
  216.    int nnl, nnr, neighbour;
  217.  
  218.    /* We run the loop from both sides until they meet */
  219.    new[0].x = old[0].x;
  220.    new[0].y = old[0].y;
  221.    new[n-1].x = old[n-1].x;
  222.    new[n-1].y = old[n-1].y;
  223.       
  224.    i = 0; j = n-1;
  225.  
  226.    while(i < j) {
  227.       
  228.       if(old[i+1].x < old[j-1].x) {
  229.      j--; do_next = j; other = i;
  230.      neighbour = i+1; nnl = i-1; nnr = i+2;
  231.      sj = 0; s = si;
  232.       }
  233.       else {
  234.      i++; do_next = i; other = j;
  235.      neighbour = j-1; nnl = j+1; nnr = j-2;
  236.      si = 0; s = sj;
  237.       }
  238. #ifdef CheapInterp
  239.       new[do_next].x = old[do_next].x;
  240.       new[do_next].y = lin_interp(new[do_next].x, old[other].x, old[other].y,
  241.                   old[neighbour].x, old[neighbour].y);
  242. #else      
  243.       new[do_next].x = old[do_next].x;
  244.       new[do_next].y = interp(new[do_next].x,
  245.                   old[nnl].x, old[nnl].y,
  246.                   old[other].x, old[other].y,
  247.                   old[neighbour].x, old[neighbour].y,
  248.                   old[nnr].x, old[nnr].y, s);
  249. #endif
  250.    }
  251.  
  252. #ifdef Debug
  253.    fputs("Interpolation debug\n", stdout);
  254.    for(i = 0; i < n; i++)
  255.       printf("%.5f %.5f   %.5f %.5f\n", new[i].x, new[i].y,
  256.          old[n-i-1].x, old[n-i-1].y);
  257. #endif   
  258. }
  259.  
  260. /*
  261.  * Very similar to get_otherside(); instead of interpolating in the other
  262.  * side, we interpolate in a different array.
  263.  *
  264.  * Caveat: when the curve turns the corner, we need to reverse the sense, so
  265.  * we fill in the points in two parts: first the upper, then the lower.
  266.  */
  267. void fill_in(point_t *interped, int n1, point_t *pts, int n2)
  268. {
  269.    int i, j=0;
  270.    
  271.    for(i = 0; i < n1; i++) {
  272.       if(pts[j+1].x > pts[j].x || interped[i+1].x > interped[i].x)
  273.      break;
  274.       while(pts[j+1].x > interped[i].x)
  275.      j++;
  276. #ifdef CheapInterp
  277.       interped[i].y = lin_interp(interped[i].x, pts[j].x, pts[j].y,
  278.                    pts[j+1].x, pts[j+1].y);
  279. #else
  280.     /*  ... *** not done yet *** ... */
  281. #endif
  282.    }
  283.    for(; i < n1; i++) {
  284.       while(pts[j+1].x < interped[i].x)
  285.      j++;
  286. #ifdef CheapInterp
  287.       interped[i].y = lin_interp(interped[i].x, pts[j].x, pts[j].y,
  288.                    pts[j+1].x, pts[j+1].y);
  289. #else
  290.     /*  ... *** not done yet *** ... */
  291. #endif
  292.    }
  293. }
  294.